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We propose a numerical method for resummation of perturbative series, which is based on the 
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the coefficients of perturbative series, and incorporates Borel resummation in a natural way. Simi- 
larly to the "worm" algorithm, the method samples open Feynman diagrams, but with an arbitrary 
number of external legs. As a test of our numerical algorithm, we study the scale dependence of the 
renormalized coupling constant in a theory of one-component scalar field with quartic interaction. 
We confirm the triviality of this theory in four and five space-time dimensions, and the instability 
of the trivial fixed point in three dimensions. 
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Introduction 

One of the most popular methods for numerical simu- 
lation of quantum field theories is the Monte-Carlo inte- 
gration over all possible field configurations, with the in- 
tegration weight being proportional to exp (—5"), where S 
is the action of the theory. However, such method might 
become not very efficient when the fields are strongly 
correlated at large distances, for example, in the vicinity 
of quantum phase transitions, or when the path-integral 
weight becomes non-positive. The latter situation is typ- 
ical, e.g., for field theories at finite chemical potential. 

Recently it has been realized that these difficulties can 
be avoided or significantly reduced if one performs the 
perturbative expansion of the theory around some "free" 
action Sq, and sums up the resulting series, instead of di- 
rectly evaluating the path integral. Such summation can 
be also performed using the Monte-Carlo procedure, and 
this alternative simulation method is usually called "Di- 
agrammatic Monte-Carlo" [1, 2]. Diagrammatic Monte- 
Carlo turned out to be very efficient for numerous prob- 
lems in statistical and condensed matter physics, and al- 
lowed to obtain many interesting results with precision 
which was unattainable for other methods. 

However, a systematic application of Diagrammatic 
Monte-Carlo to field theories which are relevant for high- 
energy physics is so far hindered mainly because expan- 
sions in powers of coupling constant, which lead to the 
conventional diagrammatic technique due to Feynman, 
typically yield only asymptotic series. Such series can- 
not be directly summed and therefore cannot be sam- 
pled by a Monte-Carlo procedure. In practical simula- 
tions, this non-convergence problem is typically avoided 
by finding a suitable strong-coupling expansion. It is 
then argued that in a finite volume this expansion can 
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be continued to the weak-coupling domain [2]. Exam- 
ples of such expansions which were used for numerical 
simulations are the strong-coupling expansion in O (N) 
and CPff lattice sigma-models and in Abelian gauge the- 
ories, and the Aizenman random current representation 
[:>] for the 0^ theory [2]. Unfortunately, the structure 
of strong-coupling expansions differs significantly for dif- 
ferent theories. Moreover, in many cases they are quite 
complicated and their structure is not sufficiently well 
understood. In particular, it has not been realized yet 
how to systematically sample strong-coupling expansion 
diagrams for non-Abelian lattice gauge theories or for 
sigma-models with SU (TV) target space. Some progress 
in this direction has been recently achieved for field the- 
ories in the large- limit [4]. 

Another way to treat the non-convergent weak- 
coupling expansions has been proposed recently in [."1]. 
The basic idea is to construct a sequence of better and 
better approximations to the original path integral, each 
of which has a convergent expansion. However, so far the 
utility of this construction was demonstrated only for the 
zero-dimensional theory. It is also not clear how to gen- 
eralize the construction of [5] to lattice field theories with 
compact variables, such as non-Abelian lattice gauge the- 
ories. One can conclude that in the context of numerical 
simulations in high-energy physics Diagrammatic Monte 
Carlo is so far a promising, but not a universal tool. 

In the present paper we propose a novel simulation 
method which, on the one hand, inherits the advanta- 
geous features of Diagrammatic Monte-Carlo, and, on 
the other hand, can be used to investigate theories with 
asymptotic weak-coupling expansions and does not re- 
quire the detailed knowledge of the structure of the 
perturbative expansion. The method is based on the 
stochastic interpretation of Schwinger-Dyson equations, 
which are typically much easier to derive than the gen- 
eral form of the coefficients of perturbative series. Such 
stochastic interpretation has been considered recently in 
[4] for large-A'^ quantum field theories, and a somewhat 
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similar approach was discussed quite a long time ago in 
[(}] (see [4] for a more detailed discussion of the methods 
used in these papers) . Similarly to the "worm algorithm" 
[ I ] , which is an essential ingredient in the Diagrammatic 
Monte-Carlo, the method samples open diagrams which 
correspond to field correlators, rather than closed dia- 
grams which correspond to the partition function of the 
theory. Another distinct feature is that the diagrams are 
sampled directly in the momentum space. Asymptotic se- 
ries can be treated within this method using the standard 
resummation tools, such as the Pade-Borel resummation 
or expansion over the basis of some special functions [7] . 
For example, in the case of factorially divergent series, 
our method incorporates the Borel transform (or, more 
generally, the Borel-Leroy transform) of the field corre- 
lators in a natural way. 

In order to illustrate the applicability of our method, 
here we consider the theory in which Schwinger-Dyson 
equations take probably the simplest nontrivial form, 
namely, the theory of a one-component scalar field with 
quartic interaction. Perturbative series in this theory are 
known to diverge factorially, and we use Pade-Borel re- 
summation to recover physical results. As a test of the 
method, we study the scale dependence of the renormal- 
ized coupling constant of the theory in different space- 
time dimensions, and confirm numerically the triviality 
of the theory in five and four space-time dimensions and 
the instability of the trivial fixed point in three dimen- 
sions 8]. For simplicity, we work only in the phase with 
unbroken Z2 symmetry. Application of our algorithm to 
the phase with broken symmetry will be discussed briefiy 
in the concluding Section, and will be studied separately 
elsewhere. We hope that the simulation strategy, which 
we illustrate here on the simplest nontrivial example, can 
be easily generalized to other quantum field theories. 

The paper is organized as follows: in Section I we intro- 
duce the basic definitions and write down the Schwinger- 
Dyson equations for the 0^ theory in momentum space. 
In Section II we provide a general description of the 
stochastic method which we use to solve these equa- 
tions. Then in Section III we formulate the simulation 
algorithm which stochastically estimates the coefficients 
of the perturbative expansion of the field correlators in 
powers of the bare coupling constant. In Section IV we 
consider a practical method to extract physical observ- 
ables from the numerical data produced by this algo- 
rithm. While Sections I - III are essential to understand 
the presented method. Section IV is more technical. In 
Section V we present and briefiy discuss some physical re- 
sults which we obtained using our algorithm. Finally, in 



Section VI we make some concluding remarks and discuss 
further extensions and generalizations of our approach. 
I. SCHWINGER-DYSON EQUATIONS FOR A 
SINGLE-COMPONENT 0^ THEORY 

We consider the theory of a one-component scalar field 
(f) (x) in Z?-dimensional Euclidean space with quartic in- 
teraction. The action of the theory is: 

5 = J d''x{l/2d^^d^(t> + ml/2(l>'' + Xo/4^^) . (1) 

The disconnected field correlators in momentum space 
are defined as follows: 



(2) 



0(p) = J d^xe'P'' (l){x) 

G{pi,...,Pn) = {(t){pi)...(f>{Pn))- 



Note that only even-order correlators arc nonzero. 
By Gc (pi, ■ ■ ■ ,Pn) we denote the corresponding con- 
nected correlators (that is, the correlators which con- 
tain only connected Feynman diagrams). Due to mo- 
mentum conservation, all correlators also contain a fac- 



tor (27r) Si ^ I . It is also convenient to define the 
\A=i ) 

correlators G (pi, . . . ,p„) from which this factor is omit- 
ted. 

We define the renormalized mass m^j and the wave 
function renormalization constant Zb. from the behavior 
of the two-point correlator at small momenta [9]: 



G'(pi,P2) 



pi+mj, + 0{pt) 



(3) 



The renormalized coupling constant A^, is related 
to the one-particle irreducible four-point correlator 
r {Pi,P2tP3tP4) at zero momenta. For the theory (1) it is 
proportional to the connected four-point correlator with 
truncated external legs: 



\r = -l/6Z^r(0, 0,0,0) = 
-1/6 Z2 (G'(0,0))"' (0,0, 0,0) 



(4) 



We note that our definition of the coupling constant dif- 
fers from the one that is used in most textbooks by a 
factor of 6. This definition is more convenient for the 
analysis of Schwinger-Dyson equations. 

Schwinger-Dyson equations for the correlators (2) in 
momentum space read: 



{2n) 



(ml+pl) G{pi,p2)^{2TrfS{p,+p2)- 
I ^^^^ ^^92 rf^ 93^ (Pi - 91 - 92 - 93)^(91,92,93,^2) 



(5) 
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{ml +pI) G {pi,p2, . ■ . ,Pn) = X! (27r)^'^(pi + Pa) G {p2, ■ ■ ■ ■.Pa-i-.Pa+i, ■ ■ ■ ,Pn) 
Ao 



(2^) 



2D 



A=2 



d qid q2d qsS {pi - qi - q2 - qs) G {qi,q2,q3,P2, ■ ■ ■ ,Pn) , 



where the arguments of the correlator in the first sum- 
mand on the r.h.s. of (6) are all the momenta except pi 
and Pa ■ 

Equations (5) and (6) were obtained by variation of 
the field correlators and the action over (p {pi ) . Similar 
equations can be obtained for any argument of field cor- 
relators pi, . . . ,pn, but the resulting system of equations 
turns out to be redundant. To obtain a complete sys- 
tem of equations, it is sufficient to consider the variation 
over (p (pi ) only. We thus arrive at a system of func- 
tional linear inhomogeneous equations for an infinite set 
of unknown functions G (pi, . . . ,p„). According to the 
general theorems of linear algebra, the solution of such 
equations is unique, if it exists. A straightforward way 
to solve equations (5), (6) is to truncate an infinite set 
of equations at some correlator order and to discretize 
the continuum momenta. In this case, however, the re- 
quired computational resources quickly grow with the 
number of correlators which arc studied. Such infinite 
systems of linear equations can be more efficiently solved 
using stochastic methods. In the next Section we pro- 
vide a general description of a specific stochastic method 
which, in our opinion, is most convenient for the solution 
of Schwinger-Dyson equations in quantum field theories. 



II. STOCHASTIC SOLUTION OF 
INHOMOGENEOUS LINEAR EQUATIONS 

We consider a system of linear inhomogeneous equa- 
tions of the following form: 

/ (.t) = ^ a (a:, y) F (x, y) f{y)+h [x) (7) 

where x is the element of some space X and ^ denotes 

summation or integration over all elements of this space. 
The coefficients A (x, y) are assumed to be positive, while 
F [x, y) can be of any sign. A [x, y) and F [x, y) arc also 
assumed to satisfy the inequalities 

Y,A{x,y) < 1 

\Fix,y)\<l (8) 

for any x^y ^ X. The source function b{x) is also as- 
sumed to be positive. Factorization of the coefficients of 
the equation (7) into A (x, y) and F{x,y) is to a large 
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extent arbitrary, and can be chosen in some optimal way 
for any particular problem. 

If the space X in (8) contains an infinite number of ele- 
ments, a deterministic approximate solution of (7) would 
require the truncation of this space to some finite number 
of elements, and it is a priori not known which elements 
can be discarded. An alternative to the deterministic 
solution is the stochastic solution, for which / (x) is pro- 
portional to the probability of occurence of the element 
X in some random process. This idea dates back to the 
works of von Neumann and Ulam (described in [lU], see 
also [11] for a more recent review). In this case, the 
solution automatically incorporates the importance sam- 
pling. Namely, those elements of space X for which / (x) 
is numerically sufficiently small are automatically trun- 
cated, since the random process cannot reach them in a 
finite number of iterations. Such methods for the solution 
of large systems of linear equations are widely used, e.g., 
in engineering and control design, for problems related 
to partitions of large graphs etc. Since the basic idea 
behind these methods is very simple, and the number of 
works on the subject is vast, we have decided to present 
here a concise formulation of the method which wc found 
most useful, without giving any specific reference to the 
literature. 

Consider the Markov process with states specified by 
one element of the space X and one real number x- We 
specify this process by the following 

Algorithm 1 At each iteration, do one of the following: 

Evolve: With probability A{x,y) change from the cur- 
rent state y to the new state x and multiply x by 
Fix,y). 

Restart: Otherwise go to the state with X = 1 '^'^'^ ^ 
random z €z X , which is distributed with probability 

m^)/(em^))- 

The "Restart" action is also used to initiate the random 
process. The first inequality in (8) ensures that the prob- 
ability of the "Evolve" action does not exceed unity, and 
the second inequality ensures that x remains bounded 
and thus have a well-defined average over the stationary 
state of the Markov process. 

Now let w {x, x) be the stationary probability distri- 
bution for this random process, that is, the probability 
to find the process in the state {x, x\ measured over 
a sufficiently large number of iterations of Algorithm 
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1. A general form of the equation governing the sta- 
tionary probability distributions of Markov processes is 
wiA) J2P{B ^ A)w{B), where P {B ^ A) is the 

B 

I 



probabiUty of transition from the state B to the state A. 
For the random process specified above such an equation 
reads: 



wix,x) = ^A {x, y) j dx 5 (x, F (x, y) x) w (y, x) 



b{x) 



6 (x, 1) ^ / dx' w (y, x') ( 1 - E ^ 



(9) 



r 



We note that Wfl = E / ^^x' w (y, x') 1 - E ^ V) 

y \ zex 

is the average probabihty of choosing the "Restart" ac- 



tion. Let us denote M ^ = wr/ \^J2 ^ (^) ) • By integrat- 
ing (9) over x one can check that 



f{x)^J\f / dxxw{x,x) 



(10) 



is the solution of the original equation (7). 

In practice, in order to find / {x) one should simulate 
the random process specified by the Algorithm 1, and 
sum up the factors x over a sufficiently large number of 
iterations separately for each x, dividing the results by a 
total number of iterations: 



Xt, 



(11) 



t=l 



where Xt and Xt are the values of x and x at t'th iteration. 
Here wr can be measured simultaneously with / (x) as 
wr = lim ni^/T, where ur is the number of "Restart" 

T-i-oo 

actions during T iterations. 



The space X should be then the space of sequences of 
momenta {pi, . . . ,p„} for any n = 2, 4, . . .. However, a 
simple analysis shows that the inequalities (8) for equa- 
tions (5) and (6) cannot be satisfied, since the number 
of summands in the first term on the r.h.s. of (6) grows 
linearly with the number n of field variables in the cor- 
relator. As will become clear from what follows, this 
growth is in fact the manifestation of the factorial di- 
vergence of the coefficients of the perturbative series. 
Such structure of Schwinger-Dyson equations for one- 
component scalar field should be contrasted with scalar 
field theories in the large- A'^ limit, where the number of 
planar diagrams grows only in geometric progression, and 
where Schwinger-Dyson equations in factorized form at 
suSiciently small coupling constant can be interpreted as 
equations for the stationary probability distributions of 
the so-called nonlinear random processes [4] . 

In order to overcome this problem, let us assume that 
the correlators (2) can be formally expanded in power 
series in the coupling constant Aq: 



G(pi, . . . ,p„) 



4-00 

m=0 



(-Ao)'" Gjn (Pl 



,Pn) (12) 



III. STOCHASTIC INTERPRETATION OF 
SCHWINGER-DYSON EQUATIONS 

At a first sight, one can try to apply the method de- 
scribed in Section II directly to equations (5) and (6). 



where the coefficients Cn,m will be specified later. 

Now we insert the expansion (12) into the Schwinger- 
Dyson equations (5) and (6) and collect the terms with 
different powers of Aq . This yields the following equations 
for the functions Gm (pi, • ■ • ,Pn): 



5,nfi (277)-° S {pi +P2) 

C2.0 {ml+pj) 

92 - q3)G,n-l (<7l, 92, 93,^2 



Gm (Pl,P2) 



+ 



C4,m — 1 



(27r)2^ iml+pj) 



C2,r. 



Jd'^qi d^g2 d'^q^Sipi-q, 



(13) 
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Gm {Pl, ■ ■ ■ ,Pn) = ^ 
Cn+2,m-l 



C„-2,m (Stt)^ 5 {pi +Pa) 



Gm {P2, ■ ■ ■ ,PA-1,PA+1, ■ ■ ■ ,Pn) + 



A=2 



{2^Y^{ml+pl) 



c„,m {rnl + p\) 

j d'^qid^q2d'^q3.5 {pi - qi - q2 - q?.)Gm-i (qi, 92, '73,P2, • ■ • ,Pn) 



(14) 



Equations (13) and (14) are also linear inhomogeneous 
equations, but on a larger functional space - the unknown 
functions G,„ (pi, . . . ,p„) now depend also on m. We can 
now try to choose the coefficients Cn,m in (12) so as to 
cast these equations in the stochastic form (7) with co- 
efficients which satisfy the inequalities (8). The space X 
now should contain sequences {pi, . . . ,p„} of n > 2 mo- 
menta in D-dimensional space and a nonnegative integer 
number m. The sum ^ should be understood as sum- 

mation over n and m and integration over n momenta: 

E^EE/'^v.-.rfV (15) 

xeX m>0n>2 

Thus, altogether the space of states of the Markov process 
that we would like to construct consists of an ordered 
sequence of momenta of arbitrary length {pi, . . . ,p„}, a 
positive integer number m and a real number x- 

Comparing the form of equations (13) and (14) with 
the general form (7), we can now construct Markov pro- 
cess which will solve these equations, as discussed in Sec- 
tion II. Again, we specify it by the following 

Algorithm 2 At each iteration, do one of the following: 

Add momenta: With probability pA = 

(2.)°(n+l)c„,„.Eo So = / 

c„+2.m ' ^ J p-'+rng ' 

\p\<A 

add a pair of momenta {p, — p} to the current 
sequence of momenta, inserting the first momenta 
at the beginning of the sequence and the second 
- between the A 'th and A + I'th elements of the 
sequence, where A is chosen at random between 
(n -|- 1) possibilities. The momentum p is dis- 
tributed within the D- dimensional sphere of radius 
A with the probability distribution ^ p2^j,i'^ ■ Do 
not change m and x- 

Create vertex: With probability pv = 
,„ ^.>ii — ^ replace the three first mo- 

menta pi, p2 and ps in the sequence by their sum 
p = P1+P2+P3- Multiply X by m^/ (ttIq + p^) and 
increase m by one. 

Restart: Otherwise restart with a sequence which con- 
tains a pair of random momenta {p, —p} with the 
probability distribution ^ p2^m'^ ; with X = 1 
and m = 0. 



Since the factor mg/ (toq -\- p^) in the "Create vertex" 
action does not exceed unity, the second inequality in 
(8) is satisfied. Let us check whether we can also satisfy 
the first inequality, that is, whether the total probability 
of "Add momenta" and "Create vertex" actions can be 
made less than one for any sequence of momenta and for 
any m. This can only be achieved if both ratios 
and — ^^^^^ — are finite for any n and m. For the first 

Cn-2,™ + l ^ 

ratio it is only possible if Cn,m grows not slower than 
(n/2)! at large n. In this case the second ratio can only 
be bounded for all n and m if c„.,„ grows as (n/2 -I- m)! 
at large n, m. We thus conclude that equations (13) and 
(14) indeed can be interpreted as the equations for the 
stationary probability distribution of some Markov pro- 
cess, if the coefficients Cn,m grow sufficiently fast. The 
factorial divergence of the perturbative series is then ab- 
sorbed in the coefficients Cn,m- It is interesting that such 
a simple analysis of Schwinger-Dyson equations reveals 
the divergence of the perturbative scries in a straightfor- 
ward way, without the need to explicitly calculate any 
diagrams! 

In Algorithm 2 we have also introduced the ultravi- 
olet cutoff A, which is necessary in order to normalize 
the probability distribution of the momenta which are 
created when the "Add momenta" action is chosen. In 
our simulations we set A = 1, so that all masses are 
measured in units of A ("lattice units"). A more self- 
consistent way to introduce such cutoff would be prob- 
ably to start with lattice theory in the coordinate space 
and assume that all momenta belong to the first Bril- 
louin zone — tt < < tt. According to the standard 
renormalization-group arguments, a particular choice of 
the cutoff prescription should result only in tiny correc- 
tions of order O (e"-'^/'^") to the physical results. There- 
fore we use the isotropic cutoff scheme, which leads to 
a much more efficient numerical algorithm. Indeed, ran- 
dom momenta with isotropic probability distribution can 
be easily generated by the standard mapping of the one- 
dimensional probability distribution w(|p|) ~ \p\'i+m'^ 
to the uniform probability distribution on the inter- 
val [0,1]. On the other hand, with lattice discretiza- 
tion in coordinate space the random momenta should 
be generated with the anisotropic probability distribu- 
tion w (pii) ^ — 'i , \^ ^ -i, — 77W. This can be done either 

by using the Metropolis algorithm or by discrctizing the 
momenta, which require much more computational re- 
sources. An interpretation of our cutoff scheme in terms 
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of Fcynman diagrams will be given below. 



is proportional to the kinematical factor 



Algorithm 2 thus stochastically samples the coefficients 
of the perturbative expansion of field correlators of the 
theory (1), which are reweightcd by the factors Cn.m- 
Each state of the Markov process defined by this algo- 
rithm corresponds therefore to some (in general, discon- 
nected) Feynman diagram with n external legs and m 
vertices. Combinatorial growth of the number of dia- 
grams is compensated by the growth of the coefficients 
Cn.m- The action "Add momenta" corresponds to adding 
a bare propagator line to the diagram, without attaching 
its legs anywhere - thus the number of external legs is 
increased by two, and the diagram order is not changed. 
The action "Create vertex" corresponds to the creation 
of one more vertex by joining three external legs of the 
diagram and attaching a new external leg to the joint. 
The number of external legs is thus decreased by two 
and the order of the diagram is increased by one. Since 
we add new momenta in pairs {p, —p}, the sum of mo- 
menta on all the external legs of any generated diagram 

n 

is always identically zero. The constraint ^ = is 

A=l 

thus automatically satisfied. 

By a slight modification of Algorithm 2, one can also 
easily trace whether the generated diagram is connected 
or disconnected. Upon the "Add momenta" action one 
should mark the two newly created legs by some unique 
label, and upon the "Create vertex" action, the labels 
of the legs which are being joined should be replaced by 
a single new label. In this way, one can directly mea- 
sure the four-point connected diagram which enters the 
definition of the renormalized coupling constant (4). 

The "Restart" action simply erases any diagram which 
was constructed before, and initiates the construction 
of a new diagram, which starts with just a single bare 
propagator. Since the maximal achievable diagram order 
is equal to the maximal number of "Create vertex" ac- 
tions between the two consecutive "Restart" actions, it is 
advantageous to minimize the probability of "Restarts" . 
The rate of "Restart" actions wr for our Algorithm 2 
thus plays the role similar to the reject probability in 
the standard Metropolis-based Monte-Carlo: for maxi- 
mal numerical efficiency it is advantageous to maximally 
reduce it. On the other hand, in order to implement im- 
portance sampling with maximal efficiency, the factor x 
in Algorithms 1 and 2 should be as close to unity as pos- 
sible. For each particular implementation of Algorithm 
1, one should find optimal balance between the rate of 
"Restart" actions and the efficiency of importance sam- 
pling by adjusting the factors F {x,y) in (7). 

In Algorithm 2 we have already chosen these factors 
to be TOq/ (toq +p^) for the "Create vertex" action, and 
unity for the "Add momenta" action. In the language 
of Feynman diagrams, such a prescription can be inter- 
preted as follows. The weight of each Feynman diagram 



Mi 



Mo 

n 



(16) 



where qi, i = 1 . . . Mj are independent momenta circu- 
lating in loops and Qj, j = 1 . . . Md can be expressed 
as some linear combinations of qi and the momenta of 
the external legs. In Algorithm 2 we in fact perform 
Monte-Carlo integration over the independent momenta 
qi, generating them randomly and independently at the 
"Add momenta" action with the probability distribution 



proportional to 1/ [qf 



Our isotropic ultraviolet 



cutoff scheme thus consists in limiting the integrations 
over independent momenta qi to \qi\ < A. The factor 
X in Algorithm 2 is then proportional to the product of 
propagators involving the dependent momenta Qj. The 
weight (16) is obtained by averaging x over random qi. 

An alternative is to perform a Metropolis-like inte- 
gration, with the probability of the "Create vertex" ac- 
tion being proportional to 1/ (p^ + TOq) and with x be- 
ing identically equal to one. However, we have found 
that such an alternative prescription increases the rate 
of "Restart" actions quite significantly, and results in a 
less optimal performance of the algorithm. Therefore, we 
consider only the first choice, implemented in Algorithm 
2. 

Let us now consider the coefficients Cn,m more closely. 
As it was already shown, Cn,m should grow as (n/2 + to)! 
at large n, to. On the other hand, if Cn,m grow too fast, 
higher-order diagrams will be strongly suppressed. Let 
us assume that c„.,„ is proportional to F (n/2 + to + a) 
times some functions which grow exponentially with to 
and n. By minimizing the rate of "Restart" actions for 
all n at a fixed to, we find the optimal value a = 1/2. 
Therefore, we choose the following form of the coefficients 



r (n/2 4- TO, + 1/2) a;-("-2)y- 



(17) 



with some x and y. The sum of the probabilities of the 
"Add momenta" and "Create vertex" actions for such a 
choice of Cm n is: 



Pa +Pv 



2 {2Trf ^qX^ (n+1) 
n -I- 2to + 1 



y 



{2n) 



2D 



.(18) 



The maximal value of the total probability is reached for 
TO = 0. Minimization with respect to x then shows that 
the total probability (18) does not exceed one only if 



y <y 



{2nf ml 



(19) 



For y = y, the total probability (18) is equal to one if 
x^ ^ x'^ ~ — . If we set y = y and a; = a;, as 

we actually did in our simulations, the probabilities of 
both "Add momenta" and "Create vertex" actions do 
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not depend on the bare mass mo or on the space-time 
dimensionahty D. Such a choice also minimizes the rate 
of "Restart" actions wn. 

Hence, the autocorrelation time of the Markov pro- 
cess specified by Algorithm 2 also does not depend on 
the physical parameters of the theory (1). In practice, it 
does not exceed several iterations. We see that our algo- 
rithm docs not suffer from the critical slowing-down in 
the sense of the usual Monte-Carlo algorithms, where au- 
tocorrelation time typically strongly increases close to the 
continuum limit. However, this absence of critical slow- 
ing down is compensated by the increase of the computa- 
tional time which is required to obtain sufficient statistics 
in the low-momentum region. Indeed, the volume of this 
region decreases as A^'^j""'^''^ as the infrared cutoff Km 
goes to zero, and hence the probability for the momenta 
Pi . . .Pn to get within this region also quickly decreases. 

It is also interesting to note that the inclusive prob- 
abilities ^ Gm (pi, • ■ • to obtain a diagram with 

rn 

n legs, irrespectively of the order to, are proportional 
to the Borel-Leroy transforms of the field correlators 
G (pi, . . . ,p„). This means that the factors (— Aq) in 
the perturbative expansion of the correlators are replaced 
by y™/r (n/2 + m + 1/2). Thus one can say that for- 
mally Algorithm 2 stochastically estimates the Borel- 
Leroy transform of field correlators for y in the range 
< y < j/. However, in order to recover physical re- 
sults, one should integrate the Borel-Leroy transform in 
the range — oo < y < 0, therefore it is difficult to use 
this statement for practical calculations. A much more 
convenient method is to analyze the dependence of the 
coefficients Gm (pi, • ■ • ,Pn) on m, as described in Section 
IV. 

For numerical experiments, we have implemented Al- 
gorithm 1 as a C program. Source code of this program 
is available at [12]. We have used the remlux random 
number generator at luxury level 2 [13]. All numerical 
data to which we refer below were obtained using this 
code. Wc have set y = y, x = x. Numerically we found 
that with such choice of parameters the average rate of 
"Restart" actions is wr = 0.282 ± 0.001. 

In order to test the performance of Algorithm 2, we 
first consider zero-dimensional scalar field theory with 
the action S {(p) = 0^/2 + A/4 0'*, for which the expansion 
coefficients in (12) can be easily obtained in an analytic 
form. On Fig. 1 we plot the ratios of the expansion coef- 
ficients Gc:n,m of the Connected n-point functions to their 
exact values. These results were calculated with 5 • 10^ 
iterations of our Algorithm 2. One can sec that up to 
approximately m = 10 numerical results differ from ex- 
act values by less than 10%, while at larger to statistical 
errors become quite significant. 

Next, on Fig. 2 we plot the probabilities of encounter- 
ing the two- and four-point diagrams of order to per one 
iteration as a function of m for the four-dimensional the- 
ory. The label "all momenta" means that the diagrams 
were added to statistics independently of the momenta of 
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FIG. 1: Ratios of the expansion coefficients Gc;n,m of the con- 
nected n-point functions for the zero-dimensional 
computed using Algorithm 2 to tlieir exact values. 
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FIG. 2: Probabilities of encountering the 2nd and 4th order 
connected diagrams in a random process specified by Algo- 
rithm 2, for all possible momenta and for momenta below 
some infrared cutoff. Solid lines are the fits of the form Aa'"^ 
with a = 0.45 



their external legs. The label p < Aj/j, where Aju = TOq 
for the two-point diagrams and A/^ = 2too for the four- 
point diagrams, means that the diagrams were addition- 



ally weighted by the factor exp 



A=l 



Pa 



2A 



(see also 



Section IV). The data was obtained for 10^ iterations of 
Algorithm 2 for = 4 and toq = 0.15. We see that when 
the diagrams are added to statistics irrespectively of the 
momenta of external legs, at large order to the proba- 
bilities of encountering both the two- and four-point di- 
agrams of order to are almost equal and decay as a~™, 
where a = 0.45 ± 0.05. This value is very close to the 
probability py = 1/2 of "Create vertex" action in Algo- 
rithm 2, since the probability of "Add momenta" action 
is suppressed at large to. Imposing the infrared cutoff 
results in a strong kinematical suppression of the proba- 
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bilities, but they still decay as a " at large m. The fits 
of the form A a~™ are plotted on Fig. 2 with solid fines. 



IV. RESUMMATION AND INFRARED LIMIT 

Algorithm 2 stochastically estimates the coefficients of 
perturbative expansion of the correlators (2), reweighted 
by the factors c„,m- In order to measure some physical 
observable one should be able to re-sum somehow the 
factorially divergent perturbative series (12). In addition, 
the zero-momentum limit of the correlators (2) should 
be taken to measure the renormalized parameters of the 
theory from (3) and (4). 

Let us first discuss the zero-momentum limit of the 
field correlators. We will first take the zero-momentum 
limit of each expansion coefficient Gm (pi, . . . ,p„) in (12), 
and then consider the resummation of the series (12). As 
discussed above, Algorithm 2 produces sequences of mo- 
menta of the external legs of Feynman diagrams, and 
Gm (0, . . . , 0) is proportional to the probability of all the 
momenta in the sequence being zero. To calculate this 
probability numerically, one should measure the proba- 
bility for all the momenta pi , . . . , p„ to belong to some 
small region near the point pi = 0, . . . ,p„ = 0, and then 
extrapolate the ratio of this probability to the volume of 
the region to the zero size of this region. Here we use soft 
infrared cutoff, and define for any function F (pi, . . . ,pn) 
of n momenta 

Fn (Air) = J d^^p^ . . . d'^p^ 

6lB,{pi,...,Pn;AlR) F{pi,...,Pn), (20) 



where 



Sir{pi, . . . ,p„;Air) 



,D/2 



v/2^ 



IR 



(«-l) 



^ exp - ^ 



A=i -^^^IR, 



(21) 



and we also assume that F (pi, . . . ,pn) contains the fac- 
tor i5 f ^ Pa] ■ This factor should be omitted when 
\A=i ) 

taking the zero-momentum limit, as in the definitions 
of the renormalized parameters (3) and (4). The zero- 
momentum limit of F (pi, . . . ,p„) corresponds then to 
the limit F„ (A//? 0). 

In order to measure the renormalized parameters 
w_R, and Z^i from the behavior of correlators (3) 
and (4), we consider the regularized zero- momentum 
limit (20) r2(A/fl,) of the quantity V {px^pi) = 
(rnj^ + Pi) G {pi , p2 ) and tune the renormalized mass m|j 
so that the deviation of r2 (A/a) from a constant value is 
minimized for A/j^ < m^j. This constant is by definition 
the wave function renormalization constant Z/j. 

The functions T2 (A/j^) are plotted for different bare 
masses and for different space-time dimensions on Fig. 3 



on the left. Horizontal solid lines are the values of Zn, 
and vertical solid lines correspond to Aju = mn. One 
can see that the limit Ajj^ ^ is indeed well-defined in 
this case, for all the values of the bare mass toq and for 
all space-time dimensions. 

In order to calculate r2 (A/^), we have to calculate first 
its expansion coefficients r2,m (A/_r), which are defined 
similarly to (12). Taking into account the definition (20), 
we measure these coefficients by summing the quantities 



I2 (m, /, Ajji) = 6n.2&iR (pi,P2; Air) (pfj 



(22) 



where I = 0, 1, separately for different diagram orders m 
over sufficiently large number of iterations of Algorithm 
2. Then 



r2,j7i (Air) 



So 



C2,0 Wr 



(^h {m,l,AiR) + m% l2 (m,0,A/fl)) , 



(23) 



where the first factor aris es due to nor malization of the 
source term in (5) and I2 {m, /, Air) denotes averag- 
ing over a sufficiently large number of iterations of the 
Markov process specified by Algorithm 2, as described in 

Similarly, in order to measure the renormalized cou- 
pling constant from (4), we should measure the one- 
particle irreducible four-point correlator F (pi,p2;?'3;?'4) 
at zero external momenta. At small external momenta 
we use (3) and (4) and write 



r(Pi,P2,P3,P4) ~ 

4 

II ^fl^ {pa + '^r)Gc{pi,P2,P3,P4)- 



(24) 



A=l 



Now one can take the regularized zero-momentum limit 
of the expansion coefficients of (24) by summing the 
quantities 

h (m, /, A/_r) 

= Sn,4SiR{pi,p2,P3,P4;AiR) £,l{pi,P2,P3,P4) (25) 

separately for connected diagrams of different orders m 
over sufficiently large number of iterations of Algorithm 
2, as in (11). Here we have defined the kinematical factors 
6 iPi,P2,P3,P4:), i = 0, ... ,4 as follows: 

Co (pi,P2,P3,P4) =pIpIpIpI 
6 {Pi,P2,P3,P4) ^pIpIpI+pIpIpI + 

I 2 2 2 I 2 2 2 
+Pl P2 P4 + Pi P2 P3 

6 {pi,P2,P3,P4) = pIpI + p\pI + pIpI + 
+pIpI+pIpI+pIp4 

6 {Pl,P2,P3,P4) ^pI+pI +pI+pI 

£.4{Pl,P2,P3,P4) =1 (26) 

Then the expansion coefficients of the regularized zero- 
momentum limit of the truncated connected four-point 
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FIG. 3: Regularized zero-momentum limits r2 (A/ij) (on the left) and r4 (A/ij) (on the right) of the two-point and four-point 
connected truncated correlators as functions of the infrared cutoff Am for Ao = 0.5 (for D = 3, we use mo = 0.15 and Ao = 0.2). 
Above: at different bare masses at D = 4, below: at mo = 0.15 in different space-time dimensions. Vertical solid lines denote 
the scale Am = m_R and the slanting and horizontal lines are linear functions which extrapolate r2 (Air) and r4 (Am) to 
Am = 0. 



correlator (24) are given by 

^4,n^ (Am) = — V mji h (m, I Km) (27) 

The resummed function r4 (A/^j) is plotted for differ- 
ent bare masses and for different space-time dimensions 
on Fig. 3 on the right (to make comparison with (4) eas- 
ier, we divide it by 6). Since the low-momentum region 
for the four-point correlator is much stronger suppressed 
kinematically then for the two-point correlator, the limit 
Kir — >■ is now not so well-defined, especially at small 
bare mass toq- In order to extrapolate r4 (A/ij) to zero, 
wc fit several (from 3 to 5) data points with smallest 
Kir by a linear function, and assume that the intercept 
of this linear function is the required limit r4 {Kir — 0). 
Slanting solid lines on the plots on the right of Fig. 3 
are these linear fits, and vertical solid lines again corre- 
spond to Kir = mR. Clearly, such extrapolation pro- 
cedure introduces quite large systematic errors into our 
measurements of the renormalizcd coupling constant (4). 



However, as we shall see below, even with such a crude 
extrapolation we get reasonable results for the scale de- 
pendence of \r. 

In order to illustrate how the soft infrared cutoff (20) 
reproduces the zero-momentum limit, on Fig. 4 we also 
show the ratios of the expansion coefficients of the reg- 
ularized zero-momentum limit of the connected n-point 
functions T2^m {Kir) and T^^rn {Kir) to the results of an- 
alytical calculations up to to = 2. Namely, we calcu- 
late analytically perturbative contributions to the two- 
point function G {p) and the connected four-point func- 
tion Gc {pi,p2,P3,P'i) up to m = 2 and multiply the ex- 
ternal legs by (to|j + p\^ , where to^ is calculated from 
our numerical data, as described below. Then we either 
set all the external momenta to zero (the corresponding 
points on the plot are shown with full symbols) or in- 
tegrate over them with the weight Sir (pi, . . . ,p„; A/i?) 
(these points are shown with empty symbols). We set 
Kir = toq /4 for the two-point function and Kir ~ toq for 
the four-point function. First we note that when both the 
analytical and the numerical results are weighted with 



10 



E 



1.2 
1.1 
1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 



■nee- 



A 



rriQ^O.IO 
mn=0.10 



^2,m. '"0^ 



mn=0.25 



^ 4, m: 



mn=0.25 



FIG. 4: Ratios of the measured expansion coeHicients 
Tn.m (A/fl) to their analytical values for different bare masses 
mo in a four-dimensional theory. For the two- and four-point 
functions we take Am — mo/4 and Aj_r = mo, respectively. 
For the data points shown with full symbols, external mo- 
menta were put to zero in analytical expressions. For data 
points shown with empty symbols, external momenta were 
integrated over with the weight Sir {pi, . . . ,pn; A/_r). 



the factor Sir{pi, . . . ,p„; Ajfj), the ratio is very close 
to one, which suggests that the loop integrals in (16) are 
indeed accurately reproduced by our Algorithm 2 for suf- 
ficiently large external momenta. On the other hand, the 
exact zero-momentum limit differs from the regularized 
result (20). As discussed above, for the four-point func- 
tion we have to take quite large IR cutoff in order to gain 
sufficient statistics, and the zero-momentum limit is then 
accessed using linear extrapolation. 

To get a deeper insight into possible problems with 
the regularized zero-momentum limit (20) of the two- 
point function, let us consider 1-particle reducible dia- 
grams with two external legs and with m insertions of 
1-particle irreducible self-energy diagrams. The contri- 
bution of such diagrams is proportional to the kinemati- 
cal factor 



(28) 



where S (p) is the self-energy. If we neglect the de- 
pendence of S (p) on p (which is indeed weak for the 
lowest-order perturbative contributions) and fix the in- 
frared cutoff in (20) to some finite A/^, at sufficiently 
large m > ttiq/Ajj^, a simple estimate shows that the 
regularized infrared limit (20) of (28) differs from the ex- 
act zero- momentum limit by a factor ^ mT'^/'^. In our 
simulations the minimal value of the infrared cutoff is 
A/fl = wo/4. According to the above arguments, this 
allows to re-sum the diagrams with up to ^ 16 inser- 
tions of self-energy diagrams. However, even for smaller 
m we expect that the contributions of the rapidly chang- 
ing kinematical factors of the form similar to (28) are 
the main source of systematic errors in our simulations. 



According to Fig. 4, they can be as large as ~ 20% for 
rriQ = 0.1 and m = 2. Presumably, these errors can be 
significantly reduced by implementing some resummation 
procedure for one-particle reducible diagrams. 

After having discussed the zero-momentum limit, let us 
describe the integral Borel-Leroy transformation which 
we use for the resummation of the series (12). For gen- 
erality, in the following we will denote by G„ any n- 
point correlator, probably multiplied by some function of 
momenta or integrated over all the momenta with some 
weight, as in (23) and (27). Correspondingly, by Gn,m 
we denote the coefficients of expansion of G„ in powers of 
the bare coupling constant Aq, reweighted by the factors 
as in (12). 

Inserting the explicit form of the factors c„.,„ from (17) 
into the series (12) and omitting the dependence on mo- 
menta, as discussed above, we obtain: 



+ 00 



Gr, 



-n+2 



^r(n/2- 



1/2) 



7n—0 



Ao 
y 



G*n,m(29) 



Using the integral representation of the gamma-function 
r (.t) = / c?f e~* i^^^ and changing the integration vari- 



able to z 



An t 



we get: 



-Tl + 2 



n + 1 +00 



Z 2 



dz exp ( — — 



Gn 



(30) 



In real simulations, the coefficients Gn.m are known 
only up to some finite maximal order. A standard re- 
summation strategy in this case is to construct the Fade 

approximant of the function G,i (z) = ^ i^^) Gn,nn 

m=0 

that is, to approximate G„ (z) by some rational function 
of z. In this case the expansion coefficients (12) are ap- 
proximated by the sum of several exponents: 



Gn 



fe 



1.1^ 
akb,. 



(31) 



where m„ is the order of the tree-level diagrams which 
contribute to the connected n-point correlator: m2 = 
and 7714 = 1- Inserting the expression (31) into the 
integral (30), we obtain: 



Gr. 



(-ir 



-n+2 



H-oo 

/ 



dz exp I — — I z 



flfc 



hkz 



(32) 



Thus each exponent in (31) corresponds to a simple pole 
of the Borel-Leroy transform. Integration over z can be 
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now performed analytically, and we obtain for the two- 
point correlator: 

G2=J2 V^^kTk (1 - x/^e'^-Erfc {^) , (33) 

k 

where = For the four-point correlator, the inte- 
gration yields: 

G4 = x^'^ ^ \f^akbk 

k 

{{rk - 1/4)' + 11/16 - e'''^Erfc (V^)) . (34) 




2 4 6 8 10 12 14 16 



m 

FIG. 5: Reweighted expansion coefficients for the regularized 
zero-momentum limit of the two- and four-point connected 
correlators T2,m (Air = mo) (23) and r4,m {Air = 2mo) (27) 
for different values of the bare mass mo in the four- 
dimensional theory. Solid lines are the fits of the data with 
the sum of exponential functions (31). Only the coefficients 
with relative error below 10% were used for fitting and are 
shown on the plot. 

The expansion coefficients F2.m {Aju ~ mo) and 
F4^,„ (A/fl = 2mo) for the four-dimensional theory, de- 
fined by (23) and (27), are plotted on Fig. 5 for different 
values of the bare mass toq- For this plot, we also omit 
the TO-independent factor (17). The renormal- 

ized mass niR used to calculate (23) and (27) corresponds 
to Ao = 0.1. Solid lines arc the fits of the form (31), with 
three exponents for the two-point correlators and with 
one exponent - for the four-point correlators. 

Due to statistical errors in the coefficients G„ the 
standard Fade approximant constructed on all the avail- 
able data points turns out to be very unstable - small 
errors in the data lead to very large deviations of the ap- 
proximant and to appearance of multiple spurious poles. 
In order to obtain more stable results, we find the opti- 
mal number of exponents and the values of ak and bk by 
fitting the expression (31) to the numerical data. 

A technical difficulty is that such fits typically result in 
badly conditioned minimization problems, if the number 



of exponents exceeds two. Fortunately, there are spe- 
cific methods which work well for this particular prob- 
lem. They are based on the singular value decomposition 
of the Hankel matrix or the so-called Page matrix [14]. 
We found that for our data the fits based on the Hankel 
matrices are optimal. Let us briefly describe this fitting 
procedure. Define the Hankel matrix Hki and the shifted 
Hankel matrix Hki as 

Hkl ~ GnM+l+m„ 
Hki = Gn,k+l+m„ + l, 

fc,/ = 0...L(m,„ax-TO„)/2j -1, (35) 

where m^ax is the maximal order to which the expansion 
coefficients are known and [. . .J is the floor function. The 
Hankel matrix can be decomposed a.s H = UTiV'^ , where 
S is the diagonal matrix with positive elements which are 
assumed to decrease from left to the right. We now form 
another matrix M = {E-'^^'^U^ HVT,-^/^) of rank N, 
with k,l — . . . N ~ 1, where N is the required number 
of exponents in (31). The eigenvalues of M are then the 
optimal values of the coefficients bk in (31) [14]. When bk 
are known, the optimal values of ak can be easily found 
by a simple linear regression. 

In our simulations we have used such maximal number 
of exponents N, for which all bk in (31) are still real and 
positive, so that all the poles of the Borel-Leroy trans- 
forms of field correlators (32) lie on the real z axis at 
z < 0. We have found that in this case the positions of 
all the poles are numerically stable. On the other hand, 
if we allow also for poles at z > or for poles off the 
real axis, their positions turn out to be numerically very 
unstable. We therefore disregard them as numerical ar- 
tifacts. There arc also theoretical arguments [15] that 
the Borel image of field correlators in bare perturbation 
theory should not have poles at positive real z. 

A disadvantage of the fitting method described above 
is that it does not take into account the statistical er- 
rors in the data when finding the optimal values of bk - 
the weights of data points do not depend on their errors. 
Therefore, after having found the optimal number of pa- 
rameters and their values in (31) from the singular value 
decomposition of the Hankel matrix, we use this values as 
an initial guess for the standard minimization-based fit- 
ting algorithm. In most cases, the optimal values which 
take into account the errors in the data turn out to be 
very close to the ones that were found from the Hankel 
matrices, and the minimization procedure quickly con- 
verges. The resulting /d.o.f is then typically of order 
of unity. We also note that we have used only the coeffi- 
cients with relative error below 10% for fitting. 

The positions of the poles of the Fade approxi- 
mants of the Borel-Leroy transforms of the functions 
r2 (A//f = Too) and r4 (A/^j = 3too) on the negative z 
axis are shown on Fig. 6 as a function of the bare mass 
toq. In order to illustrate the uncertainties in pole posi- 
tions, in Fig. 6 we show a scatter plot of poles for 10 sta- 
tistically independent data sets, each obtained with 10^ 
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FIG. 6: Positions of the poles of the Borel-Leroy transforms of the truncated connected two- and four-point correlators in the 
low-momentum region r2 (A/_r = mo) (on the left) and r4 {Air = 3mo) (on the right) on the negative z axis as a function of 
the bare mass mo, for different space-time dimensions. For better visibility, data points for different D are displaced a bit along 
the horizontal axis. 



iterations of Algorithm 2. For the two-point correlator, 
our fitting procedure reproduces three distinct poles, and 
for the four-point connected correlator, for which numer- 
ical errors are more significant, from one to two distinct 
poles can be seen. Note that when the precision of the 
numerical data is not sufficient to find two distinct poles, 
our fitting procedure yields only one pole which is situ- 
ated between the two poles which would be found if the 
precision of the data would be higher. This can be clearly 
seen for the smallest value of the bare mass toq = 0.05. 
At small bare mass the positions of the poles are also less 
stable numerically, and hence larger statistics is required 
to find them with good precision. 



V. PHYSICAL RESULTS 

In this Section we present the results of the measure- 
ments of the renormalized parameters of the theory - the 
renormalized mass m_R, the renormalized coupling con- 
stant Afl; and the wave function renormalization constant 
Zji as defined in (3) and (4). The resummation proce- 
dure and the regularized zero-momentum limit, which 
were used to obtain these results, are described in Sec- 
tion IV. All the results presented on Figs. 5-9 were 
obtained by averaging over 10 independent runs of Algo- 
rithm 2, each consisting of 10^ iterations for Z? = 3, 4 and 
4 • 10^ iterations for D = 5. For D = 5, such increase of 
the number of iterations was motivated by stronger kine- 
matical suppression of low-momentum region in higher 
space-time dimension. Averaging over independent runs 
was made in order to accurately estimate the statistical 
errors in the resummed correlators (33) and (34). Gen- 
eration of each data set with 10^ iterations took several 
hours on a single 2 GHz CPU, which is comparable to the 
computer time which was required to produce similar re- 
sults using the "worm" algorithm (several core-months 



for ^ 20 data points in different dimensions [2]). From 
Fig. 2 one can see that with such statistics from 10 to 15 
orders of perturbative expansion of field correlators can 
be analyzed in the small-momentum region (we use only 
coefficient with relative errors below 10% for resumma- 
tion, as discussed above). 
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FIG. 7: Renormalized mass niR as a function of the bare cou- 
pling constant Ao for different bare masses mo and in different 
space-time dimensions. Solid lines correspond to the one-loop 
contribution to uir. 

On Fig. 7 we illustrate the dependence of the renor- 
malized mass TTifl on the bare coupling constant Aq for 
different space-time dimensions and for different bare 
masses mp. In order to demonstrate that our resum- 
mation procedure yields the results which agree with the 
lowest orders of perturbation theory at small Aq, on Fig. 
7 we also plot the one-loop result for the renormalized 
mass, which indeed fits the data for small Aq. 

Fig. 8 shows the dependence of the wave function 
renormalization constant Zj^ on the bare coupling con- 
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FIG. 8: Wave function renormalization constant Zr as a func- 
tion of bare coupling Ao for different bare masses and in dif- 
ferent space-time dimensions. 



stant. Zfl is close to unity in the whole range of coupling 
constants and for all dimensions Z?, which agrees with the 
results of direct Monte-Carlo simulations [Ki]. Two-loop 
perturbative calculation of the self-energy shows that 
< 1, with deviation from unity not exceeding 10"'^ 
for Aq < 1.0. On the other hand, our calculations show 
that Zfi > 1 . Most likely this difference is due to system- 
atic errors in our approximation of the zero-momentum 
limit of field correlators, which were discussed above in 
detail. 




FIG. 9: Renormalized coupling constant as a function of 
renormalized mass for fixed bare coupling Aq. For D = 4 
and D = 5 we take Ao = 0.1, 0.2, 0.5, 1.0 and for D = 3, 
Ao = 0.1, 0.2. For D = 4, solid lines correspond to the result 
of integration of one-loop /3-function. For D = 5, straight 
solid lines are the linear fits of the data at small masses. For 
D = 3 and D = 5, interpolating splines are also shown to 
guide the eye. 



Finally, on Fig. 9 we illustrate the scale dependence of 
the renormalized coupling constant. Namely, we fix the 
bare coupling constant Aq (we use Aq = 0.1, 0.2, 0.5, 1.0 



for = 4, 5 and Aq = 0.1, 0.2 for D = 3) and change the 
bare mass mo . The renormalized coupling constant A r is 
then studied as a function of the renormalized mass rrifi 
in units of UV cutoff. This is, of course, equivalent to 
fixing the value of m ^ in physical units and changing the 
ultraviolet cutoff A. According to the renormalization 
group arguments [■:!, iS], when the continuum limit of the 
theory is approached and the renormalized mass mn in 
units of UV cutoff tends to zero, in space-time dimension 
_D = 4 or larger the renormalized coupling should tend 
to zero. Hence the continuum limit of the 0** theory in 
D > A dimensions is a free theory of massless scalar fields. 

Our data confirms this triviality conjecture: for both 
D = 4 and D = 5 the renormalized coupling clearly de- 
creases with the renormalized mass. For D = 4 the solid 
line on Fig. 9 is the result of integration of the one-loop 
/3-function, which implies that Xn should approach zero 
very slowly, at a logarithmic rate. Our results are con- 
sistent with such behavior within error range, although 
at small bare mass A r seems to decrease faster than log- 
arithmically. This systematic deviation from logarithmic 
scaling is probably due to large systematic errors in the 
measurement of Xr at small mg, as discussed in Section 
IV. For D — 5, at small masses the renormalized cou- 
pling goes to zero almost linearly, in agreement with the 
dimensional analysis. Linear fits of Xr on several data 
points at small renormalized mass are also shown as solid 
lines on Fig. 9. 

On the other hand, for D = 3 renormalization-group 
arguments predict that the trivial fixed point at rufj = 
0, Xn — is unstable [n]. Our data agrees with this 
statement: the renormalized coupling in this case quickly 
grows as mn goes to zero. At large renormalized mass 
ma the renormalized coupling tends to its bare value Aq 
for all bare masses and space-time dimensions. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper we have presented a novel simulation 
method which stochastically samples open Feynman di- 
agrams with probability which is proportional to their 
weight times some re-summing combinatorial factor. In 
this respect, our method is similar to the "worm" algo- 
rithm by Prokof'ev and Svistunov []_], which is often used 
in the framework of Diagrammatic Monte-Carlo. The ba- 
sic idea behind our approach is the stochastic perturba- 
tive solution of Schwinger-Dyson equations, which form 
an infinite system of linear inhomogeneous equations. 
Thus it is not necessary to know explicitly the struc- 
ture of each term in the perturbative expansion, and the 
transition probabilities are not subject to any detailed 
balance condition. In contrast to the "worm" algorithm, 
in our algorithm the number of external legs of diagrams 
is not fixed and also becomes a random variable. With 
only a minor modification of the algorithm, one can con- 
sider either disconnected or connected diagrams. 

In order to illustrate this general idea, we have applied 
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it to study the running of the renormalized coupling con- 
stant in the scalar field theory with quartic interaction. 
With our numerical algorithm we were able to obtain 
the coefficients of perturbative expansions in powers of 
the bare coupling constant up to 15tli order. Result- 
ing series were then resummcd using a specially adapted 
Pade-Borel-Leroy resummation procedure. We have con- 
firmed that the coupling constant approaches zero in the 
continuum limit in four and five space-time dimensions, 
and grows in three space-time dimensions. Performance 
of our algorithm in terms of computer time is comparable 
to the reported performance of the "worm" algorithm for 
the same theory [2]. Let us also note that the algorithm 
of [2] significantly slows down in the weak-coupling limit 
Aq — ^ 0, while the precision of our algorithm increases in 
this limit. Disadvantages of our method are the need for 
an external resummation procedure and the small signal- 
to-noise ratio in the low-momentum region. The latter 
is analogous to slowing down of the conventional Monte- 
Carlo simulations with increasing lattice volume. 

The presented algorithm clearly allows for many im- 
provements. For example, instead of expanding the cor- 
relators in powers of the bare coupling constant Aq, as in 
(12), one could try to expand them in the basis of some 
specially constructed functions of Aq (see e.g. Chapter 
16.5 of [7]). As well, one can consider other choices of 
the re-summing coefficients Cn,m in (12) than (17). For 
example, c„ „i can be proportional to the coefficients of 
the expansion of the averages ( 0" ) in powers of the cou- 
pling constant Aq in a zero-dimensional scalar field the- 
ory with the action S {(/)) = 0'^/2 + Ao 0'*/4. Preliminary 
calculations show that such choice allows to explore even 
higher orders of perturbation theory, although the resum- 
mation procedure becomes somewhat more involved. It is 
also possible to speed up the simulations by analytically 
calculating the expansion coefficients Gm {Pi, - ■ ■ ,Pn) for 
some small m. In particular, if the expansion coefficient 
Gm (pi,P2) for the two-point function is known for some 
finite TO ~ jTiQ, one can modify Algorithm 2 by start- 
ing at TO = Too and generating the first momenta in 
the sequence {^1,^2} with the probability distribution 
proportional to Gm iPi,P2)- Since to always increases in 
Algorithm 2, only diagrams of order to > too will be sam- 
pled with such a modification. With small changes, the 
method can be also applied to quantum field theories in 
the large- limit [4]. In this case the perturbative series 
are expected to diverge only exponentially, and by ana- 
lyzing the behavior of large-order expansion coefficients 
one can locate phase transitions of Gross- Witten type. 
Since the aim of the present paper is to illustrate the gen- 
eral method on the simplest possible nontrivial example, 
we do not consider here these potential improvements. 

Let us also discuss briefly the application of Algorithm 
2 to the (/)'* theory with spontaneous symmetry break- 
ing. A detailed study of this case will be presented else- 
where. Since the possibility to interpret the bare prop- 
agator as the probability distribution of the momentum 
is absolutely crucial for the formulation of Algorithm 2, 



it is not possible to set TOq < without modifying the 
algorithm. Instead, one should redefine the field vari- 
able (j) i^) ~ 4> {x) + 00 ! where 00 corresponds to one of 
the stable minima of the potential, and write Schwinger- 
Dyson equations in terms of the field (x) , which now 
has the physical mass. Note that as long as Z2 symme- 
try remains spontaneously broken, the system will not 
be able to jump to other nontrivial minimum of the po- 
tential (as it would typically happen with the standard 
Monte-Carlo in a finite volume). The reason is that our 
method in fact operates in an infinite volume limit, and 
the infrared cutoff is imposed only when collecting the 
statistics (see Section IV). 

We hope that the presented approach can be eas- 
ily generalized to other quantum field theories. It can 
be especially advantageous for those theories, for which 
the all-order perturbative expansion is difficult to con- 
struct in an explicit form or leads to asymptotic scries, 
and hence the "worm" algorithm cannot be applied in a 
straightforward way. 

Of primary interest is the extension to non-Abelian 
gauge theories with fermions, where Diagrammatic 
Monte-Carlo can potentially help to tackle the notorious 
sign problem at finite chemical potential (see, e.g., [IT]). 
In non-Abelian gauge theories Schwinger-Dyson equa- 
tions can be written in terms of gauge-invariant quan- 
tities - Wilson loops. Such formulation of Schwinger- 
Dyson equations is known as Migdal-Makeenko loop 
equations [18]. The configuration space of the random 
process which solves these equations should contain loops 
(that is, closed sequences of links) on the lattice [19], and 
the basic transformations on this space should be the 
merging and the modification of loops (for examples of 
such loop transformations see the algorithm for solving 
the Weingarten model, described in [4]). The probability 
that such random process produces some loop C should 
be then proportional to the Wilson loop W {C). The 
main problem with such approach is that straightforward 
stochastic interpretation of loop equations on the lattice 
leads to the strong-coupling expansion [4], which is not 
analytically connected with the continuum limit of the 
theory at weak coupling. Let us mention, however, that 
some attempts to explore the phase structure of QCD 
at finite chemical potential basing on the lowest-order 
strong-coupling expansion have been already reported 
in the literature, and quite interesting results were ob- 
tained [17]. Work on the apphcation of the presented 
method to loop equations in non-Abelian gauge theories 
is in progress. 
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